\documentclass[12pt,letterpaper]{article}
%\usepackage{amsmath}
\usepackage{layout, latexsym, array, enumerate, amsmath, amsthm, amssymb, amsfonts}
\usepackage[mathscr]{eucal}
\usepackage{graphicx}
\usepackage[margin=1in]{geometry}

\textwidth 7in \textheight 9.00in \oddsidemargin -0.15in
\evensidemargin -0.15in \topmargin -0.25in

\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{example}{Example}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{defn}{Definition}[section]

\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To specify the title page
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{A local-EM algorithm for spatially aggregated data with time-varying boundaries}
\author{Patrick Brown, Chun-Po Steve Fan, Jamie Stafford}
\maketitle

\section{Introduction}

Disease mapping often requires working with spatially aggregated
data, as data may be aggregated into regions to protect
confidentiality or only collected at an area level such as health
region or postal code.  These regions often vary in size, being
larger in rural areas.  Canadian six-digit postal codes are able to
locate an individual precisely in cities but rural and remote areas
have a postal codes covering large geographical extents.  Further,
estimation of disease risk must take the underlying population
distribution into account, and spatially-aggregated census data is
often used for this purpose.

Estimating disease risk for small areas or rare diseases require
using data collected over a long time period, in order to accumulate
sufficient cases to allow for accurate estimation.  As postal and
census regions change over time, this can cause problems due to the
need to combine maps with different tessellations.  Whereas a
conventional cross-sectional analysis would model risks using the
populations and case counts for census area, a longitudinal analysis
would have populations and case counts for each of the 1991, 1996,
2001, and 2006 census tracts.  Disease mapping on large, stable
areas avoids these problems, though small area analysis would
involve finer-scale census regions which are more volatile.

This paper introduces a local \textsc{em} class of algorithms for spatial data with multiple types of aggregation, {\it motivated by locally constant Poisson regression}. The algorithms extend the methods of \cite{localem} and their implementation is shown to be related to the \textsc{ems} algorithm introduced by \cite{ems}.  A local \textsc{em} algorithm is used to estimate the risk surface on a tessellation in which each of the cumulative intersections of the various maps are distinct regions. At each iteration the expected number of cases in each of the regions of the fine tessellation are computed and the risk surface re-estimated.

\subsection{Application}

The method is illustrated with an application to mapping the risk of Lupus in Toronto, Canada, from the period 1965 to 2007.  These data are collected from the lupus clinic at Toronto Western Hosptial.  Lupus may have an environmental risk factor (cite some papers) so it would be expected have a spatially structured risk pattern.

As the exact locations of the patients can be computed from the full
street address, a kernel smoother can be used to evaluate the risk
surface.  We treat this estimate as a gold standard against which we
evaluate the proposed methodology by aggregating the data to the
census dissemination area and census tract level.  Census
dissemination areas (DA) are the finest level at which population
data are released, with each region containing approximately 400
individuals and ??? DA's (perhaps) covering the greater Toronto
area.  Census tracts (CT) are larger, with ??? CT's contained in
greater Toronto.

Cancer in essex county as a (maybe) second application where exact
locations are unavailable.  partly because of confidentiality,
partly because rural addresses don't allow for precise geocoding.

\section{Methods}

The development of local likelihood and kernel methods for spatial and spatially aggregated data may be found in the literature where a sequence of papers attributed to Brillinger (1990, 91, 94)  and the methods of Silverman, Jones, Nychka and Wilson (1990), Fan and Stafford (2008) are the most relevant to efforts presented here. While Brillinger was concerned with spatially smoothing data that had been aggregated into the regions of a single map, and Silverman $et~al$ focused on the reconstruction of a single image, our methods permit multiple maps with differing tessellations to be combined. 

One view is to consider the methods developed here as an extension of the EMS algorithm to an epidemiologic setting. However the extension is motivated formally through the local-EM approach of Fan and Stafford (2008) rather than in the ad-hoc manner of Silverman $et~al$. This has the advantage of permitting 
adjustments that account for spatial variations in population size, age, sex, time, $et ~cetera$,  to arise as a natural consequence of the local likelihood construction. Points to emphasise about the proposed
methodology are that it is computationally undemanding and able to fully integrate area censoring into the estimation procedure. 

In \S 2.1 we develop a local-EM algorithm that allows multiple disease maps to be combined over time while simultaneously being spatially smoothed. Its relationship to the work of Brillinger and Silverman $et~al$ is discussed. In \S 2.2 we extend this algorithm to include an offset that permits an adjustment for other variables.

\subsection{The model}
Here disease incidence, in this case lupus within the metropolitan boundaries of Toronto, is modelled as a nonhomogeneous Poisson point process in space and time where the $i^{th}$ occurrence of the disease has location $S_i$ and occurs at time $T_i$. Initially we assume the underlying intensity $\rho$
varies spatially but not temporally, that is, $\rho(s,t)=\lambda(s)$ and interest centers on $\lambda$. The likelihood for $\lambda$ is given as
\begin{align}\label{lik}
\mathcal{L}(\lambda) &= \sum_{i} \log \rho(S_{i}, T_{i}) - \int_{\mathcal{T}} \int_{\mathcal{M}} \rho(u,v) \,
\mathrm{d}u \, \mathrm{d}v =\sum_i \log \lambda(S_i) - \sum_j |\mathcal{T}_j| \int_{\mathcal{M}} \lambda(u) \,
\mathrm{d}u 
\end{align}
where here $j$ indexes over the number of census periods in the study. For later convenience we have
written the second term as a sum over successive census periods of durations $|{\cal T}_j|$.and $\mathcal{T} = \cup_j \mathcal{T}_j$ denote the entire observation period. 
The region of study, namely metropolitan Toronto, is denoted by $\cal M$ and data is reported as counts $n_{jl}$ for regions $R_{jl}$. Here $R_{jl}$ denotes the $\ell$th census tract of the $j$th census year within the $j^{th}$ census map $M_j$ of the region $\cal M$.

Now the boundaries that define $R_{jl}$ within each map vary in time to account, for example, for the
growth of a metropolitan are as in the case of our lupus study. As a result, it is difficult to combine counts over time and to address this it is useful to construct the partition $P=\dot{\cup}J_l$ of $\cal M$ that arises from overlaying the maps $M_j$.

To estimate the intensity surface flexibly we consider the local likelihood methods of Hastie and Tibshirani (1988) and Loader (1999). Here (\ref{lik}) is replaced by
\begin{align}\label{loclik}
\mathcal{L}_s(\lambda) &= \sum_i K_h(S_i-s)\log \lambda(S_i) - \sum_j |\mathcal{T}_j| \int_{\mathcal{M}} K_h(u-s)\lambda(u) \,
\mathrm{d}u
\end{align}
and the intensity is approximated locally by 
\begin{align*}
\log \lambda(u) &=\mathcal{P}_{\bf a}(u-s)=a_0+{\bf a}_1^T({\bf u}-{\bf s})+\cdots 
\end{align*}
Here $\mathcal{P}_{\bf a}$ is a bivariate polynomial with coefficients given by ${\bf a}^T=\{a_0,{\bf a}_1,\ldots\}$ and (\ref{loclik}) becomes a likelihood for $\bf a$ at each point $s$
\begin{align}\label{loclik2}
\mathcal{L}_s({\bf a}) &= \sum_i K_h(S_i-s)\mathcal{P}_{\bf a}(S_i-s) - \sum_j |\mathcal{T}_j| \int_{\mathcal{M}} K_h(u-s)\exp\left\{ \mathcal{P}_{\bf a}(u-s) \right\}\, \mathrm{d}u
\end{align}
where the locations $S_i$ are assumed to be known. The estimate $\hat{\bf a}$ may be found by solving a set of local likelihood equations based on ${\cal L}({\bf a})$ and the estimate of the unknown intensity surface at the point $s$ is  $\hat{\lambda}(s)=e^{\hat{a}_0}$.

In our setting data is areal censored and reported as the number of occurrences $n_{jl}$ for
region $R_{jl}$. That is, the location $S$ of any particular disease occurrence is only known to be within
some region of a map. For this situation we mimic Fan and Stafford (2008) and consider an EM-type strategy where, given disease incidence is assumed to be a NHPP, the  E-step results in replacing (\ref{loclik2}) with
\begin{align}\label{loclike}
\mathcal{L}_s({\bf a}) &= \sum_{jl} n_{jl}E_\lambda \left [K_h(S-s)\mathcal{P}_{\bf a}(S-s)|S\in R_{jl}\right ] - \sum_j |\mathcal{T}_j| \int_{\mathcal{M}} K_h(u-s)\exp\left\{ \mathcal{P}_{\bf a}(u-s) \right\} \, \mathrm{d}u.
\end{align}
Here expectation is computed with respect to the conditional density
\begin{eqnarray*} \label{conden}
{{{{\lambda}}(u)}\over{\int_{R_{jl}}{{{\lambda}}(t)} \, \mathrm{d}t}}\cdot
\end{eqnarray*}
The M-step requires solving equations based on (\ref{loclike}) to get $\hat{\bf a}^r$ at the $r^{th}$ iteration. Combining the E and M steps of this local-EM algorithm 
%Now ${\cal M}=\dot{\cup} J_l$ \& $\lambda(S)=\hat{\Lambda}_{r\ell}/||J_\ell||$ for $t\in
%J_\ell$, where $\hat{\Lambda}_{r\ell}=\int_{J_\ell}
%\hat{\lambda}_r(u)\, \mathrm{d}u$. 
leads to the iteration
\begin{eqnarray}\label{localEM}
\hat{\lambda}_{r+1}(x)=\sum_{ij} \mbox{E}_{\hat{\lambda}_{r}}\left[\left.
K_h\left({S-s}\right)\right|S\in R_{ij}\right]{\Bigg /}\int_{\cal M}\tilde{K}_h(u-s)\, \mathrm{d}u
%K_h\left({S-s}\right)\right|S\in R_{ij}\right]/\Psi_h(\hat{\bf a}_r)
\end{eqnarray}
where $\tilde{K}_h(u-s)={K}_h(u-s)\exp\left\{ \mathcal{P}_{\hat{\bf a}^r}(u-s)-\hat{a}_0^r\right\}$.
%where $\Psi_h$ is
%\begin{eqnarray*}
%\Psi_h(\bf a)&=&\int_0^{\infty}K_h(u-t)\exp\left\{\sum_{j=1}^pa_j(u-t)^j \right\}\, \mathrm{d}u
%\end{eqnarray*}
%and ${\bf \hat{a}}_r$ solves a set of local likelihood equations based on (\ref{loclike}). 
Note the kernel weight for each observation is
$
\mbox{E}_{\hat{\lambda}_{r}}\left[\left.K_h\left({S-s}\right)\right|S\in R_{ij}\right]
$
and, in a strategy analogous to Fan and Stafford (2008), we simplify the computation of the kernel weight by approximating $\hat{\lambda}_r$ with a piecewise constant function $\hat{g}_r$. Here we define $\hat{g}_r(s)=\hat{\Lambda}_{r\ell}/||J_\ell||$ for $s\in J_\ell$, where $\hat{\Lambda}_{r\ell}=\int_{J_\ell}\hat{\lambda}_r(u)\, \mathrm{d}u$, so that
\begin{eqnarray*}
\mbox{E}_{\hat{\lambda}_{r}}\left[\left.K_h\left({S-s}\right)\right|S\in R_{jl}\right]
\approx \mbox{E}_{\hat{g}_{r}}\left[\left.K_h\left({S-s}\right)\right|S\in R_{jl}\right]&=&\sum_l{{\hat{\Lambda}_{rl} {\cal 
I}_{ijl}\int_{J_l}K_h\left({u-s}\right)\,
\mathrm{d}u}\over{||J_l||\sum_m \hat{\Lambda}_{rm} {\cal I}_{ijm}}}
\end{eqnarray*}
At the next iteration we are required to compute $\hat{\Lambda}_{r+1\ell}$ which leads to the simple iteration
\begin{eqnarray}\label{EMS}
\hat{\Lambda}_{r+1 s}&=&\sum_{ijl}n_{ij}{{\hat{\Lambda}_{rl} {\cal
I}_{ijl}}\over{||J_l||\sum_m \hat{\Lambda}_{rm} {\cal
I}_{ijm}}}\int_{J_s}{{\int_{J_l}K_h\left({u-t}\right)\,
\mathrm{d}u}\over{\int_{\cal M}\tilde{K}_h(u-t)\, \mathrm{d}u}}\,
\mathrm{d}t,
\end{eqnarray}
These last two expressions permit a comparison between the local-EM algorithm and methods in the literature. Note that if we only have a single map ($j=1$) then $J_l$ \& $R_{jl}$ coincide so that $I_{ijl}=0,~\forall j\neq l$. As a result the kernel weight simplifies to $\int_{J_l}K_h\left({u-t}\right)\,
\mathrm{d}u/||J_l||$, the algorithm (\ref{localEM}) iterates once and the local-EM estimator simply becomes the Nadaraya-Watson estimator advocated by Brillinger (1990, 1991, 1994) in a series of papers concerning spatial smoothing where data is aggregated to regions within a map.

Also note (\ref{EMS}) may be written as
\begin{eqnarray}\label{EMSic}
\hat{\bf \Lambda}_{r+1}={\cal M}(\hat{\bf \Lambda}_{r}) {\cal K}_h,
\end{eqnarray}
where ${\cal K}_h$ is a $k \times k$ smoothing matrix with
entries
$${\cal K}_{ls}=
{{1}\over{||J_l||}}\int_{J_s}{{\int_{J_l}K_h\left({u-t}\right)\,
\mathrm{d}u}\over{\int_{\cal M}\tilde{K}_h(u-t)\, \mathrm{d}u}}\,
\mathrm{d}t,$$ and ${\cal M}(\hat{\bf \Lambda}^{r})$ is a $k$
dimensional row vector whose $l^{th}$ entry is
$$
\sum_{ij}n_{ij}{{\hat{\Lambda}_{rl} {\cal
I}_{ijl}}\over{\sum_m \hat{\Lambda}_{rm} {\cal I}_{ijm}}}\cdot
$$
In other words (\ref{EMS}) may be written explicitly  as an EMS algorithm of the type advocated by Silverman $et~al.$ although here it is formally motivated by an EM-type strategy applied to local likelihood. Some detailed comparison of (\ref{EMS}) to Silverman $et~al$ provides further insight. The latter refers to quantities analogous to $R_{jl}$ \& $J_l$ as observation and reconstruction bins respectively. In particular, the context of Silverman $et~al$ concerns image reconstruction centered on
a single image rather than multiple maps. Nevertheless what is proposed in this paper could well be thought of as an extension of the image reconstruction techniques of Silverman $et~al.$ to an epidemiological setting. Furthermore the expression (2.2) of Silverman $et~al$ and ${\cal M}(\hat{\bf \Lambda}^{r})$ are related where, for example, their weights $p_{st}$ simplify to our indicator variables ${\cal I}_{ijl}$ because we assume the locations $S_{ij}$ have been measured without error. This observation provides an avenue for extending the local-EM toolbox to settings where data is mis-measured but this is beyond the scope of this paper.

\subsection{Offsets}

In our context we no longer assume $\rho(s,t)=\lambda(t)$ but model the time trend multiplicatively. In addition
we account for spatial and temporal variations in the underlying populations and the average incidence rate of the disease within sex and age groupings. As a result we model the intensity surface as
$$ \rho_k(x,t) = \lambda(x) {\cal O}_k(x,t). $$
where ${\cal O}_k(x,t)= \beta(t) \theta_k P_k(x,t)$. In what follws ${\cal O}_k(x,t)$ is treated as known and piecewise constant. That is, the intensity of this process is assumed to depend
on the following quantities: the population intensity $P(x,t)$, the
spatially varying relative risk $\lambda(x)$; the average incidence
rate $\theta$; and a time trend $\beta(t)$.  

To allow variations in disease rates by age and gender, consider
$P_k(x,t)$ to be a vector of population intensity by the $k$th
age-and-gender group and $\theta$ to be a vector of the corresponding disease rates.  The number of lupus incidences has a Poisson distribution, and effects are multiplicative.  That is, counts now denoted
by $n_{jkl}$, are reported by region for age and sex groups and we have
$$
N_{jk\ell} \sim \text{Poisson}\left[\rho_{jk\ell} = \int_{\mathcal{T}_j}
\int_{R_{j\ell}} \lambda(x) {\cal O}_k(x,t)dx dt
\right].
$$

Denote the indicator function $\mathcal{I}(t \in \mathcal{T}_j)$ by $\mathcal{I}_j$ and the indicator function $\mathcal{I}(t \in
\mathcal{T}_j \text{ and } x \in R_{j\ell})$ by $\mathcal{I}_{jl}$. We model discrete time effects instead of considering continous time.  Furthermore, 
assume that $P_k(x, t)$ is a piecewise constant function at each census tract during each census period.  Then we have
\[\beta(t) = \sum_j \mathcal{I}_j \beta_j \quad
\text{and} \quad P_k(x,t) = \sum_{j\ell} \mathcal{I}_{jl} P_{jk\ell}/|\!| R_{j\ell} |\!|,\] where $\beta_j$ is the $j$th period effect, and $P_{jk\ell}$ and $|\!| R_{j\ell} |\!|$ are the total population size and the area of the $\ell$th census tract at the $j$th census period. As a result we'll also have
\[\quad {\cal O}_k(x,t) = \sum_{j\ell} \mathcal{I}_{jl} {\cal O}_{jk\ell}~~ \text{where} ~~{\cal O}_{kj\ell} = \hat{\theta}_k \hat{\beta}_j
P_{kj\ell}/ |\!|R_{j\ell}|\!|\] is the offset for the $k$th age-sex
group at the region $R_{j\ell}$.
The intensity $\rho_{jk\ell}$ can then be simplified to
\begin{align}
\rho_{jk\ell} & = \int_{\mathcal{T}_j} \int_{R_{j\ell}} \lambda(x) {\cal O}_k(x,t) \, \mathrm{d}x \, \mathrm{d}t
 = \left( \theta_k \int_{\mathcal{T}_j} \beta_j  \,\mathrm{d}t \right) \int_{R_{j\ell}} \lambda(x) P_{jk}(x) \,\mathrm{d}x \notag \\
&= |\!|R_{j\ell}|\!|^{-1} \beta_j \theta_k |\mathcal{T}_j|
P_{jk\ell} \int_{R_{j\ell}} \lambda(x) \, \mathrm{d}x \notag \\
&=|\mathcal{T}_j|
{\cal O}_{jk\ell} \int_{R_{j\ell}} \lambda(x) \, \mathrm{d}x\notag \\
&=
\tilde{\cal O}_{jk\ell} \int_{R_{j\ell}} \lambda(x) \, \mathrm{d}x\label{e:
suff_stat}
\end{align}
Remarks:
\begin{enumerate}
 \item Censuses provides us with $$P_{jk\ell} = \int_{R_{j\ell}} P_{k\ell}(x) dx.$$ 
 \item The equation (\ref{e: suff_stat}) implies that $\sum_{k\ell} N_{kj\ell}$, $\sum_{j\ell} N_{kj\ell}$, and $\sum_{k} N_{kj\ell}$ are sufficient statistics for $\beta_j$, $\theta_k$, and $\int_{R_{j\ell}} \lambda(x) \, \mathrm{d}x$, respectively.
\end{enumerate}

Should the risk be spatially constant, i.e.\ no spatial variations in disease rates, we can further
simplify the integral (\ref{e: suff_stat}) to $\beta_j \theta_k
|\mathcal{T}_j| P_{jk\ell}$, and thus estimating $\theta_{k}$ and
$\beta_{j}$ using a generalized linear model (\textsc{glm}) with an
offset equal to the person-year $| \mathcal{T}_j| P_{jk\ell}$.
Next, treat the \textsc{glm} estimates of $\hat{\theta}_k$ and $\hat{\beta}_j$ as known quantities.   
To account for age and sex groupings we denote the location and time of a disease occurrence by the pair $\{S_{ik},T_{ik}\}$. Given this the likelihood function can now written as
\begin{align}
\mathcal{L}(\lambda) &= \sum_{k} \left(\sum_{i} \log
\rho_k(S_{ik}, T_{ik}) - \int_{\mathcal{T}} \int_{M} \rho_k(u, v) \,
\mathrm{d}u \, \mathrm{d}v \right) \notag\\
& =\sum_{k} \left(\sum_{i} \log
\lambda(S_{ik}) - \int_{\mathcal{T}} \int_{M} \lambda(u){\cal O}_k(u, v) \,
\mathrm{d}u \, \mathrm{d}v \right) \notag\\
& =\sum_{ik} \log
\lambda(S_{ik}) - \sum_{jkl}|\mathcal{T}_j|{\cal O}_{jk\ell} \int_{R_{j\ell}} \lambda(u)\,
\mathrm{d}u  \notag\\
& =\sum_{ik} \log
\lambda(S_{ik}) - \sum_{\ell}\tilde{\cal O}_{\ell} \int_{J_{\ell}} \lambda(u)\,
\mathrm{d}u  \label{e:llk_1}
\end{align}
where $\tilde{\cal O}_{\ell} =\sum_{jkm} \mathcal{I}(J_{\ell} \subseteq R_{jm}) \tilde{\cal O}_{jkm}$. 

The remaining details are now analogous to \S 2.1 where we initially consider modelling $\lambda$ flexibly so that (\ref{e:llk_1}) becomes
\begin{align*}
\mathcal{L}_s({\bf a}) &= \sum_{ik} K_h(S_{ik}-s)\mathcal{P}_{\bf a}(S_{ik}-s) -   \sum_{\ell}\tilde{\cal O}_{\ell} \int_{J_{\ell}}K_h(u-s)\exp\left\{ \mathcal{P}_{\bf a}(u-s) \right\}\, \mathrm{d}u.
\end{align*}
which, given the $S_{ij}$ are areal censored, is replaced by
%
\begin{align}
\mathcal{L}_s({\bf a}) & = \sum_{kj\ell} n_{kj\ell} \mathbf{E}_{\lambda} \left[ K_{h}(S - s) {\cal P}_{\bf a}(u-s) \,\Big|\, S \in R_{jl}\right] -\sum_{\ell}\tilde{\cal O}_{\ell} \int_{J_{\ell}}K_h(u-s)\exp\left\{ \mathcal{P}_{\bf a}(u-s) \right\}\, \mathrm{d}u \notag \\
%
& = \sum_{j\ell} n_{j\ell} \mathbf{E}_{\lambda} \left[ K_{h}(S - s) {\cal P}_{\bf a}(u-s) \,\Big|\, S \in R_{jl}\right] -\sum_{\ell}\tilde{\cal O}_{\ell} \int_{J_{\ell}}K_h(u-s)\exp\left\{ \mathcal{P}_{\bf a}(u-s) \right\}\, \mathrm{d}u
\end{align}
where $n_{j\ell} = \sum_k n_{kj\ell}$.  This ultimately leads to the following local-EM algorithm for estimating ${\Lambda}_l$ 
$$
\hat\Lambda_\ell^{r+1}=\sum_p \int_{J_\ell} \dfrac{\tilde{\cal O}_p \int_{J_p} K_{h}(x - s) \, \mathrm{d}x/ |\!|J_p|\!|}{\sum_n \tilde{\cal O}_n \int_{J_n} \tilde{K}_h(x - s) \, \mathrm{d}x} \, \mathrm{d}s \times %
\left(\dfrac{1}{\tilde{\cal O}_p} \sum_{jm} n_{jm} \frac{\mathcal{I}(J_p
\subseteq R_{jm}) \hat{\Lambda}_{p}^{r}}{\sum_q \mathcal{I}(J_q
\subseteq R_{jm}) \hat{\Lambda}_{q}^{r}}\right)
$$
which may again be written explicitly as an EMS algorithm
$$\hat{\bf \Lambda}_{r+1}={\cal M}(\hat{\bf \Lambda}_r){\cal K}_h.$$
Here multiplication and division by the quantity $\tilde{\cal O}_p$ ensures ${\cal M}(\hat{\bf \lambda}_r)$ is a step in an EM algorithm and also that ${\cal K}_h$ has rowsums equal to 1 {\em add comment that this matters for convergence}.
Finally, upon the convergence an \textsc{ems} estimate for the relative risk $\lambda(x)$
is
\begin{equation}
\hat\lambda^{\ast}(x) = \sum_p \dfrac{\tilde{\cal O}_p \int_{C_p} K_{h}(x - s) \,
\mathrm{d}x/ |\!|C_p|\!|}{\sum_n \tilde{\cal O}_n \int_{C_n} \tilde{K}_h(x - s) \,
\mathrm{d}x} \times %
\left( \dfrac{1}{\tilde{\cal O}_p} \sum_{j\ell} \frac{\mathcal{I}(C_p \subseteq
R_{j\ell}) N_{j\ell} \hat{\Lambda}_{p}^{\ast}}{\sum_q
\mathcal{I}(C_q \subseteq R_{j\ell}) \hat{\Lambda}_{q}^{\ast}}
\right)
\end{equation}

Note that the matrix expression (\ref{e: ems}) is recognized as an
\textsc{ems} algorithm, in which $\mathcal{M}$ is an \textsc{em}
mapping with an extra smoothing step represented by $\mathcal{S}_h$.

% This local likelihood is then, which leads to a local \textsc{em} algorithm.

%likelihood for $\lambda$
%$\begin{align*}
%$\log L(\lambda) =& \sum_k \sum_{i; X_i = k} \log \rho_k(s_i, t_i) - \int \rho_k(s, t)\\
%= & \sum_k \sum_j \sum_{\ell=1}^{M_j} \sum_{i; X_i = k, T_i,S_i \in R_{j\ell}}
%\end{align*}
%
%only depends on
%
%\[
%O_{j\ell} = \sum_k \theta_k P_{jk\ell} \beta_j (t_j - t_{j-1} )  / |R_{j\ell}|
%\]
%and the $S_i$ and $T_i$, not the $X_i$ .
%
%
%$S_i$ are poisson point process, intensity = ?
%
%$\lambda(x) = \lambda$, constant in space, MLE
%
%local likelihood, uniform kernel and general kernel
%
%
%Area censoring results in the $S_{it}$ being unobserved and we are only provided with the number of cases $Y_{tj}$ in region $R_{tj}$ during the $t$th period. Writing $P_{tj}$ to be the vector of populations for each age and sex group in region $j$ and period $t$, we assume the population is evenly distributed within each region. Although this assumption is not ideal, census regions and $\theta$ to be the population-average disease rates for each age-sex group, we



\subsection{Simulation Results}

A simulation study of 100 samples is conducted to assess the performance of the proposed local EM relative risk estimator in the simplified scenario that assumes only one age-and-sex group and no temporal variation in the spatial risk surface. Each simulated sample consists of two maps that represent two observations in the same geographical area at different time points. Although the outer boundaries of these two map are identical, the subregions of these two maps are of different shapes. More specifically, there are five vertically stacked horizontal rectangles in the first map and five horizontally juxtaposed vertical rectangles in the second map.  All subregions have the same area of 5 squared-units. Overlaying these two maps forms 25 unit-squared cells on which kernel weights are computed.  A Gaussian kernel and a sequence of 201 equal-distant bandwidths, which ranges from 0 to 2, are chosen to construct the kernel weight matrix as shown in (\ref{e: ems}).

Once the boundaries of the subregions are defined, spatial Poisson point process data are simulated in the following steps.  First, the population locations are simulated using a Poisson process whose intensity function is uniform over each subregion.  The population intensities per squared-unit for the five subregions in each map are set to be 18, 28, 38, 28, and 18.  The true relative risk $\lambda(x_i)$ is set to be the product of rescaled two gamma density functions with the shape and scale parameters of 1.5 and 0.5, respectively.  This relative risk function attains the maximum of 100 at $(0.25, 0.25)$.  Next, for each simulated observation $x_i$, a case label is randomly generated with the probability of $\lambda(x_i)/100$.  Finally, the population and case data are both aggregated over the subregion where they are located and reported in the form of regional counts.

We compute two relative risk estimates for each simulated sample and for each bandwidth.  One relative risk estimate is the proposed EMS estimate, and the other is a kernel version of the EM estimate smoothed by placing the estimated massess at the centres of the squared cells.  In addition, the kernel intensity estimate for exact locations is also computed.  Here, we choose the mean integrated squared errors (MISE) to compare these relative risk estimates.  For each risk estimate and each bandwidth, the MISE is estimated by averaging the integrated squared errors of the 100 simulated samples and plotted in Figure (\ref{f: sim_mise}).

As shown in Figure (\ref{f: sim_mise}), the proposed EMS risk estimate attain a lower MISE than the smoothed EM risk estimate, and the lowest MISE is closer to the one when observations are exact.  Moreover, the bandwidth that the local EM achieves the lowest MISE is considerably smaller than that of the smoothed EM one.

\begin{center}
\vspace{-1.5cm}
\begin{figure}
% \includegraphics[scale=.6, angle=0]{mise_simulation_v2_gaussian}
\caption{The proposed local EM intensity estimate achieves the lowest
overall MISE with a small bandwidth of 0.195, comparing to the smoothed EM
estimate by placing expected increments at the centres of pixels.}
\label{f: sim_mise}
\end{figure}
\end{center}



\section{Application}


\subsection{Cross Validation}

Leave-one-map-out cross validation is used to choose an optimal bandwidth size.  Since excluding a map alters the offsets $O_{n} = \sum_{j\ell} \mathcal{I}(C_{n} \subseteq R_{j\ell}) O_{j\ell}$ in the $(p, m)$th entriy of $\mathcal{S}_h$, the cross validation requires the re-calculation of the smoothing matrix each time a map is left out, which leads to four smoothing matrices of 206,957 by 206,957 for each bandwidth value. To facilitate this computationally intensive task, we limit the construction of the smoothing matrix to 14 equal-distant bandwidth values, ranging from 150 to 2100.  Moreover, the geographical areas that these four maps represent are different because the Greater Toronto Area is expanding during this 40-year observation period.  We restrict the estimation of prediction error, which is the squared difference between the observed and predicted values,  over the census tracts that are common to all four maps, and there are 2,689 common census tract.  Specifically, $$ \mbox{PE}(h) = \frac{1}{4} \sum_{j\ell'} \left( N_{j\ell'} - O_{j\ell'}\hat{\Lambda}_{j\ell'}(h) \right)^2,$$ where $\ell'$ is the index of the common regions.  The bandwidth that gives the smallest predictioin error, is chosen to be the optimal, and the optimal bandwdith value 1350.

\section{Discussion}


\section{Computational Details}
Although the EMS algorithm is easy to implement, iterative procedures become computational demanding and time consuming when
the number of pixels is large.  To reduce this computational burden, we takes the following two steps to expedite the necessary computational process.  First, we choose a radially symmetrical kernel with compact support to construct the smoothing matrix $\mathcal{S}_h$, which results in a sparse matrix, i.e.\ most matrix entries equal to zero, when the bandwidth is small relative to the pixel and map sizes. In this case, we can apply algebraic procedures for sparse matrices to speed up the EMS iterations with minimal computational resources.

Next, we need to numerically evaluate the quadruple integrals that
appear in the smoothing matrix $\mathcal{S}_h$.  A brute-force
method is to apply the Gaussian quadrature rule to the four
dimensional function $$\int_{C_m} \dfrac{\int_{C_p} K_{h}(x - s) \,
\mathrm{d}x/ |\!|C_p|\!|}{\sum_n O_n \int_{C_n} K_h(x - s) \,
\mathrm{d}x} \, \mathrm{d}s;$$ however, this method becomes
computationally infeasible when the number of pixels is large.
The computational burden of these quadruple integrals can
be greatly reduced by deriving an analytical solution in a closed form
to the inner double integrals and then apply the Gaussian quadrature
rule to approximate the outer double integrals. For these two reasons, we choose the bivariate biweight function as
the kernel function
$$K(u, v) = \begin{cases}
         \dfrac{3}{\pi} \left(1-(u^2 + v^2)\right)^2 & \mbox{where $u^2 + v^2 \le 1 $}\\
        0   & \mbox{otherwise}
        \end{cases}
$$ in the construction of the smoothing matrix $\mathcal{S}_h$. 

For any fixed $s$, we then apply Green's Theorem to simplify the double integral $\int_{C_p} K_h(x-s)\,\mathrm{d}x$ to a line integral, thus calculating this double integral exactly. For instance, let $h=1$, $u=x_1 - s_1$ and $v=x_2 - s_2$. Denote the intersection of the unit circle and the pixel over which the integral is computed by $C$ and the boundary of $C$ in a counter clockwise orientation by $\partial C$. Then the inner double integral is given by
\begin{align}
 & \int_{C} K(u, v)\,\mathrm{d}u\,\mathrm{d}v %= \frac{3}{\pi} \int_{C} \left(1-(u^2 + v^2)\right)^2 \,\mathrm{d}u \,\mathrm{d}v \notag \\
= \frac{3}{\pi} \int_{\partial C} f(u, v) \,\mathrm{d}u + g(u, v)\,\mathrm{d}v, \label{e:line_int}
\end{align}
where $f(u, v) = -(v^5/5 - 2v^3/3 + u^2v^3/3+v/2)$ and $g(u, v) =
u^5/5 -2u^3/3 +u^3v^2/3 + u/2$. In terms of polar coordinate with
$r=1$, the equation (\ref{e:line_int}) can be expressed as
\begin{align}
 \frac{3}{\pi} \int_{\partial C} (5 \theta - \sin 4\theta)/30 \, \mathrm{d}\theta
\end{align}
%
When $\partial C = \bigcup_i \partial C_i$, where $\partial C_i$ is
either the boundary of the pixel or the arc of the unit circle,
$\int_{\partial C} f(u, v) \,\mathrm{d}u + g(u, v)\,\mathrm{d}v =
\sum_i \int_{\partial C_i} f(u, v) \,\mathrm{d}u + g(u,
v)\,\mathrm{d}v$ and
 $$ \int_{\partial C_i} f(u, v) \,\mathrm{d}u + g(u, v)\,\mathrm{d}v =
    \begin{cases}
    \frac{3}{\pi} \int_{u_R}^{u_L} f(u, v) \,\mathrm{d}u & \mbox{if $\partial C_i \equiv (u_R, v) \rightarrow (u_L, v)$} \\
    \frac{3}{\pi} \int_{v_U}^{v_L} g(u, v)  \,\mathrm{d}v & \mbox{if $\partial C_i \equiv (u, v_U) \rightarrow (u, v_L)$} \\
    \frac{1}{10\pi} \int_{\theta_L}^{\theta_U} 5 \theta - \sin 4\theta \, \mathrm{d}\theta  & \mbox{if $\partial C_i \equiv \theta_L \rightarrow \theta_U$.}
    \end{cases}
$$
%
Once the inner double integral is evaluated for any fixed $s$, the
outer double integral can be then numerically evaluated using the
Gaussian quadrature rule, i.e.\
$$
\int_{C_m} \dfrac{\int_{C_p} K_{h}(x - s) \, \mathrm{d}x/ |\!|C_p|\!|}{\sum_n O_n \int_{C_n} K_h(x - s) \, \mathrm{d}x} \, \mathrm{d}s \approxeq \sum_i w_i \dfrac{\int_{C_p} K_{h}(x - s_{mi}) \, \mathrm{d}x/ |\!|C_p|\!|}{\sum_n O_n \int_{C_n} K_h(x - s_{mi}) \, \mathrm{d}x},
$$ where $s_{mi}$ is the $i$th quadrature point in $C_m$ with quadrature weight $w_i$.



\end{document}
